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Abstract 

We review the long term project of the ALPHA collaboration to compute in QCD the running coupling constant 
and quark masses at high energy scales in terms of low energy hadronic quantities. The adapted techniques required 
to numerically carry out the required multiscale non-perturbative calculation with our special emphasis on the control 
of systematic errors are summarized. The complete results in the two dynamical flavor approximation are reviewed 
and an outlook is given on the ongoing three flavor extension of the programme with improved target precision. 
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1. Introduction 

Quantum Chromo Dynamics (QCD) is the renormal- 
izable quantum field theory containing gluon and quark 
fields that interact in a unique way dictated by SU(3) 
gauge invariance. It may be seen as arising from the 
standard model of elementary particles in a limit where 
all other fields, including their interactions with quarks 
and gluons, are stripped away. Strong interactions and 
confinement are the characteristics of this sector which 
hence calls for non-perturbative evaluations and is in the 
focus of lattice formulations and simulations. 

We here consider QCD with a free number of Nf color 
triplets (flavors) of quark species. In Nature we see the 
case Ni - 6 with the flavors up, down, strange, charm, 
bottom and top in order of ascending mass. The species 
beyond light up and down quarks come with charac¬ 
teristic scales of the order of 0.1 GeV, 1 GeV, 4 GeV, 
175 GeV. Therefore it makes sense to consider effec¬ 
tive theories with Vf < 6 to describe physics with char¬ 
acteristic energies significantly below the scales of the 
dropped degrees of freedom. They then enter only in¬ 
directly into the determination of the free parameters of 
the effective theory. In lattice simulations the modelling 


of the precise flavor content is technically very demand¬ 
ing. Therefore a lot of studies are found and will also 
be discussed here that refer to Af = 2 and Nf -h where 
the latter number is the minimum to allow for real ap¬ 
plications as an effective theory H] |2l [3 [H |5l. In any 
case generalized QCD has Af -i- 1 free parameters given 
by one quark mass per species and in addition the gauge 
coupling. 

The two light species are in most studies, including 
those described here, approximated to be degenerate. 
The value zero for some or even all Af quark masses 
is theoretically nice as it enhances the chiral symme¬ 
try of the model and is thus stabilized under renormal¬ 
ization. The renormalization of the coupling can be 
defined in this massless limit and we then speak of a 
massless renormalization scheme. Such schemes are 
technically convenient in nontrivial perturbative as well 
as non-perturbative calculations. The renormalization 
of the coupling can be left unchanged as quark masses 
are ‘turned on later’. To define a renormalized coupling 
constant in a massless scheme an additional scale p en¬ 
ters via the renormalization conditions. The resulting 
scale dependent ‘running’ coupling gip) obeys a Callan- 
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Symanzik equation 

^ Pigiff)), ( 1 ) 

dll 

in which the function y6 is determined by the theory once 
a particular coupling definition has been adopted. A 
negative yS-function corresponds to asymptotic freedom. 
The free integration constant that arises in solving this 
differential equation can be taken as the free parameter 
that the bare coupling has been ‘traded for’ in the pro¬ 
cess of renormalization. It may be fixed by specifying 
g for a specific yu value in GeV. Alternatively one may 
convert the Callan-Symanzik equation into the equiva¬ 
lent integral statement that 
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is independent of yu. In this equation bo, b\ are the lead¬ 
ing and scheme independent coefficients in the asymp¬ 
totic expansion 


n>0 
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For asymptotically large n it is sufficient to evaluate Q 
with the perturbative series for (3 truncated beyond some 
n > 1. Therefore, in a perturbative context, A is associ¬ 
ated with the behavior of g{p) for p ^ oo. 


2 . Hadronic renormalization scheme and finite size 

scaling 

In lattice simulations also non-perturbative quantities 
associated with scales of order one GeV and below can 
be computed in principle. Examples are the masses of 
light hadrons and matrix elements involving their one- 
particle states, decay constants like f„,fK for example 
iiiii. This opens up the possibility to also match such 
quantities directly to experiment and in this way deter¬ 
mine the free parameters of QCD which can then be de¬ 
termined with an in principle arbitrary precisiorQ This 


* This refers to pure QCD. Other interactions are still neglected. 


is not true if perturbation theory at any finite energy is 
involved, since with an asymptotic expansion - even 
if very high orders were available - an uncertainty re¬ 
mains. This effect is expected to be small at the Z-mass, 
but the situation is much more delicate for example for 
determinations of otj in the r-mass region. 

As a conceptually simple example of a hadronic 
scheme one could imagine to use as input parameters 
the mass of the proton and in addition the masses of 
Nf types of stable mesons that are sensitive to the re¬ 
spective quark masses. In practice one of course has 
Nf + I dimensionless parameters at ones disposal in the 
lattice theory of which Nf may be determined by di¬ 
aling the correct ratios of meson to proton mass. The 
remaining degree of freedom allows to tune the lattice 
theory to its critical point where the continuum limit is 
reached. Due to asymptotic freedom in QCD this is ac¬ 
complished by sending the bare coupling to zero. In this 
limit, all dimensionfull quantities emerge in the form of 
well-defined multiples of appropriate powers of the pro¬ 
ton mass which we thus employ to set the scale for all 
observables. Equivalently we may say that all that is 
computed from theories including the lattice and com¬ 
pared with experiment are dimensionless ratios of ob¬ 
servables. The above scheme selects a minimal set of 
independent mass ratios and, with these tuned, all other 
ratios must ‘fall in place’. We try to be very explicit 
on this seemingly trivial issue, as sometimes confusion 
seems to arise here nevertheless. 

In the previous paragraphs we have described a rather 
idealized situation. Eor various technical reasons we 
will not use this precise hadronic scheme, and in ad¬ 
dition several sources of in practice unavoidable sys¬ 
tematic errors have to be taken into account in lattice 
computations. 

A lattice that is simulated on a computer necessarily 
has a finite number of sites and thus finitely many de¬ 
grees of freedom. This implies a finite extent L and a 
finite spacing or resolution a such that one has (Lja)* 
sites. In large present day simulations Lja ~ 100 is 
achieved. If we refer to mhad as some hadronic mass 
scale, then a/Whad > 0 represents a distortion of the 
physics by an unphysical UV cutoff effect. Details de¬ 
pend on the chosen lattice discretization, but in prac¬ 
tice and employing Symanzik’s theory of cutoff effects 
El nmol, we expect these effects to diminish asymp¬ 
totically at a rate proportional to (awihad)^- We need 
to verify that we have reached this asymptotic behav¬ 
ior to estimate the prefactor by multiple simulations in 
which the resolution (and nothing else) is varied. This 
whole procedure is called continuum extrapolation and, 
of course, leaves behind a contribution in the final error 
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budget. 

In the same way, unless finite size effects are de¬ 
liberately looked for (see below), also a finite product 
i-'Whad < oo is an unwanted IR cutoff effect to be extrap¬ 
olated away or at least bounded. The situation here is 
however more benign, as theory ifTTll implies that the in¬ 
finite volume limit is reached at an exponential rate pro¬ 
portional to exp(-m;rL), where the pion with mass m,j 
enters as the lightest degree of freedom. In todays large 
volume lattice simulations we are largely restricted to 
cutoffs L < 6 fm and < 5 GeV, although these 
extreme values cannot yet quite be realized simultane¬ 
ously and compromises, depending on the physics stud¬ 
ied, have to be made. 

One would clearly like to confront results extracted 
by matching perturbation theory to experiment at high 
enough energy with those of a hadronic scheme. This 
amounts to nothing less than establishing QCD as one 
theory at all these length scales. An obvious strategy 
is to compute the renormalized parameters of a pertur¬ 
bative scheme from ‘within’ a hadronic scheme. Fo¬ 
cussing on the coupling constant this would require to 
compute beside a hadron mass mhad some scale depen¬ 
dent observable Oip) that possesses a perturbative ex¬ 
pansion 

0(41) = a(n) + pia^Qj) + pio^ip) + ... ( 6 ) 

where the coupling a refers to some perturbative 
scheme like MS. This evaluation has to proceed at a 
high enough scale ///whad = p » 1 to get sufficient per¬ 
turbative precision in extracting a at energy p x Whad- 
With the help of Q this information may be converted 
to a value for A/mhad- As is well known, the A param¬ 
eter of any other scheme follows now by relating the 
corresponding couplings at one loop order. 

If we imagine to perform such a calculation naively 
on the lattice, we have to cope with a multiscale problem 
where we have to satisfy the string of inequalities 

a « « L. (7) 

Given the practical constraints on L/a it is clear that 
such a direct approach will require severe compromises. 
An overview over various approaches including the di¬ 
rect one is found in Q. 

Due to these difficulties strategies have been devised 
to alleviate the problem by circumventing one of the re¬ 
quired large scale ratios. One idea is to try to tolerate the 
scales a ' and p to lie in the same range and thus per¬ 
form perturbation theory at the scale of the cutoff. Such 
a calculation, as well as references to earlier versions, 
is discussed in HU. A valid criticism in our opinion 


is, that it thus becomes hard to disentangle UV cutoff 
effects from limitations of (truncated) perturbation the¬ 
ory. To control lattice artefacts it seems necessary to 
be able to vary the lattice spacing over some range with 
physical scales mhad,p held fixed. 

Our finite size strategy which is in the focus of the re¬ 
mainder of this article may be seen as identifying p and 
L * over a major part of the calculation. We will exploit 
the fact that finite size effects are universal predictions 
of quantum field theory. The Casimir force in QED lfT3l 
is such a finite size effect which has been experimentally 
confirmed as a subtle manifestation of vacuum fluctua¬ 
tions. We extend this concept to more abstract cases in¬ 
volving periodic or Dirichlet type boundary conditions 
on finite systems. They are not realized in the labora¬ 
tory, but universality here means that there again are 
unique predictions for effects depending on the size L, 
independent of how the field theory is regularized and 
the UV cutoff limit in the finite box is taken. 

We now sketch such a computation and come back to 
details in the later sections. The key quantity required 
in our approach is an L-dependent finite size observable 
g{L) that can function as a non-perturbatively defined 
coupling constant. It must exist for arbitrary L and pos¬ 
sess a manageable perturbative expansion that is appli¬ 
cable at small size, for L ' much larger than hadronic 
scales. With a suitable normalization g{L) will be re¬ 
lated to other couplings like for example 

f{L)^t-(p) + c(pLrg^-ip) + 0(gy) ( 8 ) 

where for good perturbative accuracy one will take 
pL - 0(1). Beside perturbation theory giL) must be 
easily computable numerically and a reasonable sig¬ 
nal to noise ratio for its estimator is another practical 
requirement. Appropriate such couplings are known 
for Schrodinger functional boundary conditions, see the 
later sections for details. 

With g{L) at hand we now first consider a simula¬ 
tion in a large, i. e. effectively infinite, volume where 
we tune the bare parameters to achieve a « 

L. A hadronic scale m\^, or rather the dimensionless 
combination a/Whad is the output. Then, for the same 
bare parameters we diminish Lja to the point where 
Lmhad = (T/a)(a»*had) = P becomes a fixed number of 
order unity (the previous value amhad is used). At this 
size g{L) will not be small but its value can be computed 
by simulation. By repeating these steps for several small 
a/Whad and extrapolating to the continuum we thus de¬ 
rive a numerical value g{L) at the scale L - pmf^^ in 
the non-perturbative regime. This is a universal number 
^(Tmax = P^hld) continuum where L^ax can be 

cited in fermi once mhad has been related to experiment. 
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What remains to be done now is to evolve g(L) from 
the now known starting point at L^ax to other L small 
enough to make contact with couplings like by per¬ 
turbation theory. The Callan Symanzik equation is a dif¬ 
ferential description of such an evolution guided by the 
corresponding yS-function. For a non-perturbative evo¬ 
lution the change of scale by a finite factor two seems 
more natural. This leads to the definition of a finite step- 
size counterpart of a y6-function, the step scaling func¬ 
tion IflTlI 

o-(m) = f(2L)\gi(L)^u- (9) 

Note that cr does not refer to the lattice at this point 
but is a universal continuum quantity, different however 
for differently defined couplings. Before we outline the 
computation of the function cr we discuss its applica¬ 
tion. We use it to build a sequence {m,, i - 0,1,... ,n) 
based on the recursion m,- = cr(ui+i) and started at 
Mo = g^(Tmax)- This immediately implies that 

( 10 ) 

holds and that for sufficiently many steps this value is 
arbitrarily deep in the perturbative regime. Here in ad¬ 
dition we should also start to see an evolution that co¬ 
incides with the one produced by the Callan Symanzik 
equation with a perturbative approximation for the ap¬ 
propriate jS function. In figure a schematic view is 
given for the example of g realized by the Schrodinger 
functional finite volume coupling that will be detailed 
below. 



S(2.u,l/6) 



Figure 2: Illustration of the computation of the continuum step scaling 
function from finer and finer lattices. Note that 'L{2,u,alL) in the 
illustration corresponds to our lattice step scaling function Ii(w, a/L). 

the bare parameters such that g^ - u results and that the 
quark masses vanish (massless scheme). Then, keeping 
the bare parameters fixed, we double the size to 2L/a 
and determine a value of the lattice step scaling function 
S(m, fl/L) = g^(2L). If this procedure is repeated for a 
sequence of Lja we finally extrapolate to 


= O(lfm) : HS SF(p = 1/L„„,) 

i 

SF(p = 2/L„„,) 

i 


cr(u) - lim l,(u,alL) (11) 

a/L-^0 

which again is a piece of universal continuum physics. 
We need to implement this at a sufficiently dense set of u 
values to have a sufficiently precise control over cr(u) in 
the range required for the evolution ( [T0| ). The procedure 
is indicated in figurej^for one value of u. 


I 

SF(p = 2"/Lmax) 
PT: i 

DIS, jet-physics, at s = <— Aqcd 

Figure 1: The strategy for a non-perturbative computation of short 
distance parameters. SF refers to the Schrodinger functional renor¬ 
malization scheme and HS to a hadronic scheme. 

It remains to outline the computation of cr. To this 
end we pick a resolution Lja and some value u and tune 


3. Lattice discretization 

3.1. The action 

The gluon vector potential A^{x) in Euclidean contin¬ 
uum QCD has its values in the Lie algebra SU(3) and 
enters as a parallel transporter over infinitesimal dis¬ 
tances into the covariant derivative - d^j + Af^. The 
resulting gauge covariant curvature F^y - {D^,Dy'\, also 
in the Lie algebra, is the building block of the gauge in¬ 
variant action density that is integrated over in the non- 
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( 12 ) 


Following Wilson Ea we put strong emphasis on 
gauge invariance and construct the discretization on a 
hypercubic lattice in a way that is manifestly compati¬ 
ble with this structure. The smallest separation on the 
lattice are finite links (x,j(i) starting from any site x into 
the ju-direction and ending at its nearest neighbor x + afi 
(a is the lattice spacing and /) a unit vector). The associ¬ 
ated gauge transporters U(x,fj.) are SU(3) group-valued 
and represent the fundamental gluon field on the lat¬ 
tice. In the path integral they are integrated over with 
the invariant Haar measure independently on each link. 
Curvature on the lattice shows up by the two-step trans¬ 
porters U(x,fj.)U(x+ap.,v) and U(x,v)U(x+av,fj.) being 
different or, equivalently, by the plaquette field 


P(x,ju,v) = 
U(x,fi)Uix - 


(13) 

afi, v)[U{x, v)U(x + av,ii)Y^ 


differing from unity. Thus a direct transcription of (12 1 
onto the lattice is given by the Wilson plaquette action 


5w ^ Retr[l -P{x,ii,v)]. 


(14) 


x,p<v 


If we build U(x,fj.) out of slowly varying infinitesimal 
continuum A^(x) we find that S w approaches 5cont with 


= 6 /gQ, i. e. (|14| is a classically valid discretization. 


In the same rdference Wilson has also introduced 
what today are called Wilson lattice fermions, which 
will be the discretization of the quark Dirac field em¬ 
ployed here. Independent quark fields ij/{x) are intro¬ 
duced on the sites of the lattice. They carry a four¬ 
valued Dirac index, the SU(3) color index and an A^f 
valued flavor index which we all suppress in our nota¬ 
tion. The covariant forward and backward derivatives of 
quarks are defined by 

(D^i//)(x) = -[U(x,Y)Hx + afi) - il/{x)] (15) 


(D*i/^)(x) = -U (x- a(i,n)il/{x - ajl)] (16) 

It is well known that the most obvious discretized Dirac 
operator leads to spurious particle poles, the so-called 
fermion doubling problem. Wilson’s remedy for this 
problem is the addition of a term proportional to a dis¬ 
cretized form of the covariant Laplacian, the Wilson 


^ In our convention and are antihermitean. 


term, which attributes a mass of order a * to the doubler 
modes but is otherwise suppressed by an extra power of 
a. The total Wilson Dirac operator is now given by 




(17) 


2 ^ ^ 2 
The total lattice QCD action is thus given by 

S w;w[17, i/t, 1 ^] - S w[17] - ^ iA(Dw + (18) 

;c 

where Mq is the bare mass matrix with the mass param¬ 
eters for the various flavors on its diagonal. All these in¬ 
gredients may now be Anally assembled to write down 
the lattice QCD partition function 


f 


Z= DUDi/^Di/^expi-S W;w[f7, l/r, l/r]) 


(19) 


which beside the group integrations requires additional 
integrations over the Grassmann valued quark fields. 
For numerical simulations the latter (Gaussian) integrals 
are carried out and produce the fermion determinant 


f 


Z= Dt/exp{-5w[t/])det(Dw+ Mo). 


( 20 ) 


The matrix Dw has t/ fields in its matrix elements and 
the determinant represents a complicated nonlocal ef¬ 
fective action in U that has to be taken into account 
for sampling t/ with the weight given by this integrand. 
The standard hybrid Monte Carlo algorithms are able 
to cope with this, but nevertheless here is the source of 
the enhanced complexity of simulations once dynami¬ 
cal quark degrees of freedom are included. Some more 
details on our algorithmic implementation of the QCD 
path integral will be given in section]^ 

Another important point to mention is that the Wilson 
fermion regularization breaks chiral symmetry which 
emerges only in the continuum limit. Due to the Wilson 
term, Dw does not anticommute with 75 . As a conse¬ 
quence the masses on the diagonal of Mq undergo addi¬ 
tive renormalization and the physical zero mass condi¬ 
tion to set up a massless scheme has to be enforced as a 
nontrivial renormalization condition - usually some chi¬ 
ral Ward identity - which will force aMg to approach a 
go dependent nontrivial critical value amdgo)- Note that 
in perturbation theory one finds arndgo) - cigo + O(.?o^ 
which amounts to a linearly diverging bare mass param¬ 
eter. 

3.2. Symanzik improved action 

A further consequence of missing chiral symmetry 
in the regularized theory is the appearance of cutoff ef¬ 
fects that are linear in the lattice spacing a (multiplied 
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by powers of logarithms). In Symanzik’s approach the 
structure of cutoff effects can be studied by describing 
the lattice theory including leading cutoff effects by an 
effective theory in the continuum. This requires addi¬ 
tional terms in the action of the latter beyond the combi¬ 
nation of renormalizable terms of dimension up to four 
that we have discretized before. In asymptotically free 
theories the dimension of terms can be used to organize 
the additional contributions to the action; extra terms 
of dimension five represent lattice artefacts that vanish 
linearly in a, dimension six those that are quadratic and 
so on. The second essential criterion is to only admit 
terms that are invariant under the (reduced) symmetry 
that is still present in the lattice theory. Here chiral 
symmetry is missing for Wilson fermions. The inter¬ 
play between dimension and symmetry leads to a finite 
number of additional couplings that allows to match all 
cutoff effects up to a given order in a. This is the stan¬ 
dard situation in effective theories, with a rapid prolif¬ 
eration of the number of terms as the order is increased. 
In addition, all that one can hope for is an asymptotic 
expansion that is relevant close to the continuum limit 
where combinations like amhad are already small. In ad¬ 
dition to the enlarged action also in renormalized ob¬ 
servables extra mixings with higher dimensional terms 
of the same symmetry have to be taken into account 
to achieve a complete representation of cutoff effects in 
correlation functions. The whole concept may be seen 
as an extension of the renormalization programme that 
normally just focuses on divergences (log a and possibly 
inverse powers) to small positive powers of a. 

The dimension five terms relevant for the linear order 
in a have been classified and listed in ifl^ for mass de¬ 
generate quarks. In this case there are five operators of 
which three can be absorbed into modifications of the 
bare coupling and the bare mass. The remaining two 
operators are the Pauli term 

Ol = (21) 

and 

O 2 = -H (22) 

Note that within the Symanzik effective theory we sys¬ 
tematically expand in a. Therefore the dimension five 
terms in the action with an explicit factor a are expanded 
down from the exponent and appear as operator inser¬ 
tions. 

The description and ultimately elimination of a class 
of cutoff effects becomes much more manageable if we 
restrict ourselves to what is called on-shell improvement 
m. The restriction refers to correlations with a finite 


number of local operators which all reside at physical 
separations from each other, finite multiples of n\l^ for 
example. Then we can deform the integration variables 
of the path integral at all points without inserted opera¬ 
tors to derive the so-called equations of motion as oper¬ 
ator identities. They can be used to transform between 
the higher dimension terms classified before and to re¬ 
duce their number while still matching cutoff effects in 
on-shell correlations. For the case at hand the equation 
of motion is just the Euclidean Dirac equation 

{y^D^ + Mo)ili = 0, - Mo) = 0 (23) 

where only the renormalizable terms are considered as 
we want to only modify the order five terms and ne¬ 
glect yet higher dimensional terms. This equation is 
usually employed to eliminate O 2 and then the only 
new bulk term needed to describe 0(a) on-shell cut¬ 
off effects in the effective theory is the Pauli term Oi. 
One more complication has to be mentioned. The inser¬ 
tions generated by the a-expanded improved action ap¬ 
pear integrated over Euclidean space-time and thus are 
not strictly separated from the observables as the use 
of the equations of motion would require. A more de¬ 
tailed analysis confirms however that these violations of 
the equations of motion due to overlapping insertions 
(‘contact terms’) can be compensated in the observable 
improvement terms that were mentioned before. 

The Symanzik description of cutoff effects can be 
used in a next step to eliminate these contributions by 
what is called Symanzik improvement. To that end a 
discretized version of the extra term(s) is added to the 
original action with coefficients that are tuned such that 
the corresponding couplings in the effective action van¬ 
ish. This may then be seen as an alternative discretiza¬ 
tion without the leading artefacts that have been sys¬ 
tematically canceled in this way. Eor the operator Oi 
above this is the so called clover term first proposed in 
ifTTl . We then have the Sheikholeslami-Wohlert (SW) 
improved lattice action 

S W;SW = S W;W + ^ ^(x)-^CrftyFpyix) lA(x) (24) 

X 

where the name clover is owed to the lattice represen¬ 
tative of the field strength F^y{x) from four co-planar 
plaquettes (‘leafs’) 

Fpvix) = -^{Q^iyix) - Qyffx)], (25) 

Qpvix) - P(v,/j, v) H-three more (26) 

rotated 1x1 loops opened at x. 
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see figure An improvement condition (generalized 
renormalization condition) that has to hold in the con¬ 
tinuum theory has to be enforced to determine csw(^o)- 


Xv 


Figure 3: The clover leaf representation F^v(-r) of Fyjv(x). 

The number of terms and coefficients required at 
leading order is small enough that we can and shall im¬ 
plement complete on-shell 0(a) improvement. This is 
not practicable any more at the next order a^. If one im¬ 
plements however only some part of the improvement 
terms with some prescription for the coefficients one 
still obtains a legal variant discretization different from 
the one without extra terms. Numerical experience sug¬ 
gests that the addition of a rectangle term to the gluon 
plaquette action with a strength suggested by improving 
at tree level of perturbation theory leads to a variant ac¬ 
tion with better properties than the pure plaquette form 
although artefacts are O(a^) in both cases. This action, 
called tree level improved Liischer-Weisz action, gener¬ 
alizes ( [T4| ) and reads 

1 

5lw = 

i=0 CeSi 

Here So is the set of all different (unoriented) plaquette 
(1x1) loops on the lattice and PiC) a parallel transporter 
around it. Hence, for co = l,ci = 0 this action would 
coincide with 5w- The second term involves the set .Si 
of all different planar 1x2 loops (rectangles) and for the 
tree level improved Liischer-Weisz action the weights 
assume the values 

« = 5 , C, = - 1 . (28) 

3.3. Improved currents and renormalization 

Quark currents are observables of primary impor¬ 
tance in QCD. In particular the isovector axial current 
formed from the two light quarks 



(29) 

and the pseudo-scalar density 


P^(x) = iftys-P'ift(x) 

(30) 


enter into the discussion of chiral symmetry. They are 
here given first as bare currents in terms of bare fields 
at the same lattice site and Pauli matrices operate on 
the up and down quarks. As discussed in the previous 
subsection for 0{a) improvement these dimension three 
operators can mix with dimension four terms of the right 
symmetry. It turns out that there is no such term for P", 
but the improved axial current can mix with the gradient 
of P" and is hence given by 

(Ai);(x) = A“(x) + flCA (31) 

Here 5^ = (5^ -H 5p/2 is the symmetrized lattice deriva¬ 
tive and CA(go) is an improvement coefficient that has to 
be fixed by another improvement condition. It will turn 
out that its perturbative expansion starts at ©(g^). 

In USl the nontrivial interplay between the use of a 
massless renormalization scheme and improvement is 
discussed in some detail. In such a scheme all renor¬ 
malization conditions are formulated at a normalization 
scale p. For nonzero, but for simplicity degenerate, 
quark masses mo the relation between bare and renor¬ 
malized coupling must be taken as 

lo = ^o(l + bgamq). (32) 

Here bg(go) is an improvement constant that eliminates 
0 (a) effects at nonzero m^ which in turn is the sub¬ 
tracted quark mass 

mg = mo - mdgo) (33) 

such that niq - 0 implies a vanishing physical mass. 
The term with bg reflects a dimension five term in the 
Symanzik effective action proportional to mtr{Fj)y). It 
is only with this term (and the correct bg) that in the 
process of expressing physical observables in terms of 
gR not only divergences but also linear lattice artefacts 
are eliminated. 

In a similar way the usual multiplicative mass renor¬ 
malization must be replaced by 

mR = ihgZ^igl, ap), fhg ^ mg{\ + b,„amg). (34) 

Quite similar formulas follow for the current renormal¬ 
izations 

(Ar)“ = ZA(l+^Aflm«)(Ai);, (35) 

(Pr)“ = Zp{\+bvamg)P^. (36) 

If in the continuum limit the chiral symmetry group 
SU(Af) X SU(Af) is recovered up to finite mass effects, 
then we expect the PCAC relation 

5^(Ar); = 2 mR(PR)‘= 


(37) 
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space 

(LxLxL box with periodic b.c.) 

Figure 4: The Schrodinger functional geometry. 

to emerge as an operator relation that can be inserted 
into matrix elements. It is very advantageous to focus 
on this relation to actually define the renormalized mass 
in the improved theory, once the improved renormalized 
currents have been introduced already. The finite size 
Schrodinger functional scheme that we introduce next 
is a very convenient setting to do so. 

Non-perturbative techniques to determine the im¬ 
provement coefficients were developed in the quenched 
approximation lfT8lfT9ll20l . and later applied for the two 
and three flavor theories in 11211122ll2^l24l . Some coef¬ 
ficients such as by remain unknown non-perturbatively 
but can be taken from one-loop perturbation theory 

mM- 

4. The Schrodinger functional 

For the Schrodinger functional (SF) ETlI - and later 
the associated renormalization scheme - we consider a 
finite portion of Euclidean space time with spatial ex¬ 
tent L and temporal size T as depicted in figure]^ In 
the spatial directions fj. - k - 1,2,3 we impose periodic 
boundary conditions x=x±Lk while Dirichlet bound¬ 
ary conditions fix certain field components at xq - Q,T 
to externally given values. The SF can be defined in 
the continuum and is in fact studied by dimensionally 
regularized 1-loop perturbation theory in ETll . As we 
here want to mainly discuss non-perturbative computa¬ 
tions we prefer to immediately start on the lattice, where 
some features of the SF even become simpler to discuss. 
This implies, of course, that both Lja and Tja must be 
integer. 

4.1. Gauge sector 

In the standard form of the SF the Dirichlet condi¬ 
tions for the gluon field U{x,/i) are 

= exp(aC^), (38) 

t/(-T,k)U„=r = exp(aQ, (39) 


where are constant Abelian vector potentials in 

the form of diagonal traceless imaginary matrices. Note 
that temporal links U{x, 0) exist as integration variables 
for .To - 0,a,... ,T-a and are not subject to any bound¬ 
ary conditions. 

A possible interpretation of this Euclidean path inte¬ 
gral with boundaries is in the Hamiltonian or transfer 
matrix formalism. It goes with a Schrodinger represen¬ 
tation of states in Hilbert space as wave function(al)s on 
(three dimensional) spatial configurations U{x,k). We 
denote a state concentrated on a fixed field configura¬ 
tion U{x,k) = exp(aCk) by a ket |C) (like |x)-states in 
quantum mechanics). Then the SE partition function X 
is equal to the matrix element 

Z(C'; C) = <C'|e^™P|C), (40) 

where e""® is the transfer matrix and P is a projector to 
gauge invariant states ll28l . 

Erom the Euclidean point of view in the SE setup the 
time direction is distinguished from the others and there 
is no translation invariance in the time-direction. There¬ 
fore in this case we generalize the Wilson action to 

S Wsf = w(C)Retr[l-P(C)]. (41) 

CcSq 

Here the novelty is the plaquette dependent weight w 
for which the transfer matrix formalism suggests to take 
w{C) - 1 for all plaquettes except the purely spatial 
ones on the boundary where w(C) - 1/2. This is so 
because it is natural to symmetrically distribute these 
contributions to the two adjacent transfer matrix factors 
(as for the potential V{x) in quantum mechanics). 

In the Symanzik effective action additional terms rep¬ 
resenting cutoff effects are possible due to the SE geom¬ 
etry. In the continuum they are given by three dimen¬ 
sional integrals over boundary planes with the dimen¬ 
sion four densities triFokFok), ^(FkiFki) as integrands. 
These contributions are associated with 0(a) boundary 
cutoff effects. Eor Symanzik improvement these terms 
are discretized and included in the lattice action with 
adjustable coefficients. Eor 5wsf this is incorporated 
by two different nontrivial weights for space-time and 
for space-space plaquettes at the boundary. With our 
Abelian boundary fields the latter type does not con¬ 
tribute and we set 

w(C)-c,(go) = l+c''Vo + --- (42) 

for Ok plaquettes touching the boundary and w(C) = 1 
for all others. In principle c,(go) has to be fixed by 
yet another improvement condition. In practice these 
terms can at present only be set to perturbative values 
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as already indicated above. See ll^ for a discussion of 
the relevance of this approximation. For the tree-level 
Liischer Weisz action similarly modified weights are re¬ 
quired for both plaquettes and rectangles close to the 
boundary which are discussed in IIMII or more recently 
in ll3TII . In the following we always assume that these 
Oia) improvement terms are included in lattice actions 
for the SF. 

To not introduce further scales the boundary fields 
C, C' are taken as multiples of L * and an often used 
standard choice is 

Ct = ^diag(77-7r/3,-77/2,-77/2+ ;r/3) (43) 

C[ = ^diag (-77 - ; 7 r, 77/2 H- 7 r/ 3 , ? 7 / 2 -H 2;r/3). 

The dimensionless parameter 77 allows to vary the 
boundary values and is set to zero after taking deriva¬ 
tives with respect to it. In Il27l it is shown that these par¬ 
ticular boundary values (for not too small Lja, T/a) lead 
to a stable minimum of S wsf■ This minimum is unique 
up to gauge transformations which (at the boundaries) 
we restrict to the subgroup that preserves the boundary 
values. A representative - exp(aB^) for this 

gauge orbit of minima is 

Bo = 0, = [xoC, + (T - xo)Ck]/T (44) 

which linearly interpolates between the boundaries. 

We are now in a position to introduce the renormal¬ 
ized SF coupling. We start from the effective action or 
free energy 

F[B] =-lnZ(C';C). (45) 


In perturbation theory a saddle point expansion around 
B requires the usual Fadeev-Popov gauge fixing and 
yields a regular expansion 

r[B] = 5wsf[B] + ri[B] + glT 2 [B] + ..., (46) 


where we note that the classical or tree level term 
‘5wsf[B] in (41 1 is proportional - y6/6. The deh- 
nition of the coupling associated with the scale L finally 
reads 


-Ins- 2 ^‘^Wsf /577 


^sf(T) = ^0 


dTIdrj 


(47) 


,7=0 


In dTIdrj we differentiate the logarithm of a partition 
function with respect to a parameter entering into the 
(boundary terms of the) action. It immediately leads to 
an expectation value that is independent of unphysical 
factors in the path integral measure and can be estimated 
as a mean value of a well defined observable when 
field configurations are sampled with the probability 
exp(-5wsf)- The normalization factor g^dSwafldr] can 
easily be computed in closed form for ([44|. 


4.2. Quark sector 

The extension of the SF to fermions has been first 
presented in 1321 . 

We generalize the periodic boundary conditions in 
space to a periodicity up to a phas^for the quark helds 

i/'(x + Lk) - e'®*i/7(x), ij/{x + Lk) — e^'®‘i/'(x), (48) 


which leaves all bilinear densities in the action and else¬ 
where strictly periodic. The angles Ok allows us to vary 
the finite size kinematics in useful ways. 

As for the Dirichlet boundary conditions in time, it 
turns out that only half of the components of the in¬ 
dependent Grassmann fields i/r and i/r have to be fixed. 
Formally this is a consequence of the hrst order nature 
of the Dirac equation and the correspondingly modified 
boundary value problem as was already noted in con¬ 
nection with the bag model ll33l . Alternatively, and 
more laboriously, one may argue on the basis of the 
transfer matrix for Wilson fermions 01. The result for 
the SF in any case are Dirichlet conditions 

= P(x), P-iAlxo=r =p'(x) (49) 

and 


*pP-\xa=Q = P(X), = P'(X) 


(50) 


with projectors 

P± = ^(i±ro)- 


(51) 


The spatial helds p,p',p,p' are formal Grassmann val¬ 
ued sources which, after possible differentiations will 


always be set to zero. The partition function (45 1 now 


depends on all boundary helds Z, - Z(C', p', p'; C,p,p) 
and is given by the hnite volume path integral 

X- \ Dt/Di/7D(/7exp[-5w;Wsf(17,(A,lA)] (52) 


which as an example we have written with the plaque- 
tte action and plain Wilson fermions and the (here sup¬ 
pressed) boundary helds enter into the action. 

As mentioned before we want to eliminate 0(a) arte¬ 
facts throughout. We know already that this requires the 
clover term to be added to the bulk quark action. But, 
as for the gluons, the presence of boundaries allows for 
new terms in the Symanzik action that have to be can¬ 
celed by corresponding improvement terms. A discus¬ 
sion of the relevant dimension four densities integrated 


^ It has become customary to call this twisted boundary conditions, 
although it should not be confused with t’Hooft type twisted boundary 
conditions refening to planes of a torus. 
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shows that differentiation with respect to p(x) for ex¬ 
ample leads to single insertions of the dynamical quark 
fields with suitable parallel transporters. They transform 
contragrediently to p(x), i. e. like p(x). Hence the var¬ 
ious ■ can be contracted to form boundary currents 
like for example 

O" (57) 

u,v 


Figure 5: Correlation functions in the Schrodinger functional. On the 
left we show boundary-to-bulk correlation functions such as /p and on 
the right the boundary-to-boundary correlation function f\. 


over the xq - Q,T time slices is found in ifThl together 
with their reduction due to the on-shell conditions. The 
upshot is a rather simple modification of the quark ac¬ 
tion by a contribution dl xq - a given by 

(q - \)a^ (53) 

X 

and another such term with the same coefficient at the 
other boundary. Again c, has to be tuned and possesses 
a perturbative expansion 

Q(go) = 1 + c''Vo + • ■ ■ ■ (54) 

Additional terms involving spatial boundary terms do 
not contribute for constant sources p,..., to which 
our applications will be restricted, or can be absorbed 
into rescaling the sources which undergo multiplicative 
renormalization anyway. If in the following we refer 
to the Sheikholeslami Wohlert improved quark action 
S swsf for the SF, we assume the c, terms to be included, 
too. 


4.3. Boundary quark operators 

We introduce ‘functional’ differentiation operators 
for the boundary sources 


r(x) 


-3 d 
“ dp(x) ’ 
-3 -9 

5p'(x)’ 




-3 


dp(x) 


(55) 


r(x) = -a 


dp'(x)' 


If they contribute to an observable O expectation values 
are meant in the sense 


= (56) 

DUDif/Dij/O exp[-5 (U, if/, if/)} I 

j p'=p'=p=p=0 

where S is any of our lattice actions and it contains 
the boundary values. Inspection of any of these actions 



and the corresponding primed operator at the other 
boundary. They allow to form the SF standard corre¬ 
lation functions 

Mxo)^-^-{P%x)0‘‘} (58) 

and the boundary to boundary constant 

/i - (59) 

which are illustrated in figure]^ We may use these cor¬ 
relation functions to define a convenient normalization 
condition for Zp in ( [36] l by postulating 

Zp = const. s/TilfpiT 12) (60) 


where multiplicative renormalization factors of the 
boundary fields cancel. Some choice has to be adopted 
for the kinematical parameters: aspect ratio T/L, the 
sources, angles Ok and for xq/T. The constant is then 
chosen such that Zp = 1 holds at tree level. As em¬ 
phasized before we would like to establish the SF as a 
massless renormalization scheme. We therefore want to 
tune the bare masses mo to their critical value. At least 
in perturbation theory this is possible as the SF supplies 
an infrared regulator by providing a mass gap of order 
L-\ 

A similar standard matrix element involving the axial 
current is defined by 

fA{xQ)^-^-{Al{x)(r) (61) 

which allows us to define a bare improved PCAC mass 


dofA(xo) + CAad^do/pixo) 
2/p(xo) 


(62) 


In the next step this leads to the renormalized running 
mass 


m(L) — m 


Za(1 -bAam,f) 

Zp(l - bpaniq) 


(63) 


If all improvement and renormalization factors assume 
their correct values m is well defined up to O(a^) un¬ 
certainties. This is visible as small variations under 
changes of kinematic parameters (xo for example) on 
which the lattice but not the continuum results depend. 
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5. Gradient flow 

The finite volume SF coupling defined in the previous 
section has been the key tool for the ALPHA collabo¬ 
ration’s effort to relate low and high energies in QCD 
for many years. Recently the technique of the gradient 
flow ll35l has given access to a large class of alterna¬ 
tive observables that may serve in the same function. It 
has turned out in the meantime that it is a good strat¬ 
egy to keep SF boundary conditions, typically with zero 
boundary fields, to define massless finite size schemes 
based on gradient flow (GF) couplings 0^ . They will 
be found superior in practise to the SF coupling at in¬ 
termediate energy scales, but not at the perturbative end 
of the scale evolution. In addition, perturbation theory 
is more developed and also easier for the traditional SF 
coupling where the three-loop term of the yS-function as 
well as some two loop results for improvement coeffi¬ 
cients like c, are available (for S w;SWsf) IITH . Therefore 
a good way forward is to use a particular optimized GF 
coupling for the lower energy part of the scale evolu¬ 
tion and, at an intermediate energy, to match it to the 
SF coupling which is then run up to perturbatively high 
energies. This promises the best over-all precision and 
is thus a good reason to also exploit gradient flow cou¬ 
plings and review them here. 

5.7. Gradient flow in the continuum theory 

In the continuum Yang Mills theory we follow lf35l 
and grow a one parameter (f) family of gauge potentials 
B^{t, x) starting in an arbitrary four dimensional poten¬ 
tial A/jix) which is taken as an initial value 


the dimension of a squared length. The five dimensional 
action is constructed with additional Lagrange multipli¬ 
ers in such a way that the flow equation appears as a 
constraint in the functional integral. Similar techniques 
have been known for some time in the related field of 
stochastic quantization If39ll . 

A particularly striking new finite observable is the lo¬ 
cal action density 

E(t)^- ^tr(G^v(f, x)G^,y{t, x)), (66) 

where G^y(t, x) is the SU(3) curvature tensor of the po¬ 
tential Bft(t,x), which has a finite expectation value at 
positive flow time. In llTSl the leading order perturba¬ 
tion theory result 




has been derived. To leading order the length is the 
radius over which the initial values A^(x) are smoothed 
by solving (65 i and is thus a natural scale for the MS 
coupling to expand in. This shows that up to a com¬ 
putable normalization factor the combination fl{E{t)} 
can be regarded as a renormalized coupling constant 
which - in contrast to - is however defined also 
beyond perturbation theory. In section we have ex¬ 
plained the advantages for our purposes in defining cou¬ 
pling constants running with the finite size. This re¬ 
mains true for GF couplings. We want however to stay 
with couplings running with only a single scale and 
therefore tie together f, L by setting 


^t^cL 


( 68 ) 


Bfjit = 0, x) = Aij(x). (64) 

The trajectory is defined by imposing the first order flow 
equation 

dBflt,x) _ ^6 Sym[B] 

dt 6B^{t,x) 

where 5 ym is the usual Euclidean Yang Mills action 
evaluated for potentials at any t value, which will be 
referred to as flow ‘time’. The action is lowered with 
growing t as we move along the steepest descent direc¬ 
tion of the action. Thus along the trajectory the initial 
field gets smoothed. 

It is most remarkable that correlation functions at fi¬ 
nite t are finite without additional renormalizations be¬ 
yond those of the four dimensional theory. In ll38l this 
has been shown to all orders of perturbation theory by 
performing an analysis based on Feynman rules of a five 
dimensional theory where t is an extra coordinate with 


with a fixed proportionality constant c of order unity. Its 
value and the details of the finite size boundary condi¬ 
tions together form part of the coupling definition and 
each choice represents a different scheme. These differ¬ 
ent variants are in finite relations with each other and a 
particular choice will be determined by practical consid¬ 
erations. In QCD with Nf quark species the flow ( |65] l is 
unchanged with only S ym appearing. Also the leading 
order result ( [67| is Nf independent. 

We remark that the result ( [67l i indicates that GF cou¬ 
plings are not technically simple in perturbation the¬ 
ory. The normalization factor in defining a renormalized 
coupling constant via the bare coupling is usually given 
by a trivial tree level result. For GF couplings a dif¬ 
fusion process corresponding to the free flow equation 
has to be solved to normalize E{t) and this calculation 
is already comparable to the evaluation of one-loop dia¬ 
grams. The next correction would resemble a two loop 
calculation etc. 
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5.2. Gradient flow on a finite lattice 

In principle the transcription to the lattice is straight 
forward, but many choices have to be made that differ in 
their lattice artefacts. First of all, a lattice discretization 
for the action in the four dimensional path integral has to 
be chosen as before, with which configurations 
are generated. In addition, the action whose gradient is 
taken in the flow equation must be discretized, and this 
need not necessarily be the same discretization. Finally, 
a lattice substitute for the action density must be picked. 

On the lattice the gluon field is Lie group rather than 


algebra valued. This requires standard changes to (651. 
Instead of B^i{t,x) we now introduce a family of group 
valued link field V{t, x,/j.) by 


V{t - 0, x,/i) = U(x,ii) 


(69) 


and 

2dV{t,X,)l) 1 2 a CTT/n t^r^s 

a -- V(t,x,fi) ^-god,.c,fjS[V]. (70) 

Here S is any lattice gluon action (e. g. 5 w or 5 lw)- The 
Lie algebra valued left derivative gradient of a scalar 
function on the group is given by (suppressing here the 
link index x, fl) 


df[U] = J]r^/[exp(5nt/] 


(71) 


with r" being a basis of group generators. Such a lattice 
version of the gradient flow based on U(x,fi) fields has 
been first written down in llLSll . 

In ll40l a careful assessment of lattice artefacts of gra¬ 
dient flow observables has been reported. We do not re¬ 
view any details which proceed via an analysis of the 
five dimensional Symanzik effective theory. The re¬ 
markably simple result is that the major part of all O(a^) 
coming from the flow equation is cancelled if we in¬ 
voke the tree level improved Liischer Weisz action for 
the gradient flow with just one additional term and re¬ 
place ([70|) by 


2 dV(t, X,Ll) 


V(t,x,fi)-^ = 


(72) 


-gl(\ + , 


where the covariant derivatives defined in ( [T5] l, ( [Thl l 
have to be taken with V(t,x,fl) here. The discretized 
energy density to form E has to be read off from the 
density in 5 lw <|27]l. It has to be remembered however 


that the usual 0((?) effects of the four dimensional the¬ 
ory are unchanged, and it is just the additional source 


from the flow that one tries to control here. More nu¬ 
merical experiments with the rather recent proposal ( [72) i 
are required. 

To evaluate flow observables in simulations some 
solver for first order equations in time has to be em¬ 
ployed. This will require also a discretization of flow 
time, but this step-size error can be kept at a negligible 
level. In appendix C of 051 a third order Runge Kutta 
integrator on Lie groups has been proposed in a fully ex¬ 
plicit form that is easy to implement. It is reported that 
a step-size e - 0.01 in tja^ is accompanied by errors 
of 10“^ in link variables. A more sophisticated variant 
in 0 ^ automatizes the choice of e by combining the 
previous Runge Kutta integrator with adaptive step-size 
control. 

To define a finite size scheme based on a GF coupling 
given by the expectation value of {E{t)} boundary con¬ 
ditions have to be specified. The first attempt to do so 
was made for simple periodic boundary conditions in 
ED. It has been known for a long time that Yang Mills 
theory on a small torus is complicated due to the pres¬ 
ence of non-Gaussian fluctuation modes (see refs, in 
El). This leads to a small coupling expansion of torus 
based finite size couplings that is non-analytic in g^. As 
discussed in section]^ the SF in contrast has a unique 
minimum with purely Gaussian fluctuations up to exact 
gauge modes that have to be fixed in the usual way and 
then no non-analyticity arises. Moreover, zero field SF 
Dirichlet boundary conditions provide an infrared reg¬ 
ulator and hence the additional bonus is that the quark 
masses can be set to zero and we continue to have a 
massless scheme as with the SF coupling before. 

To achieve the exact normalization on the lattice 
N^^fl{E(t,xo)) - g^ + 0(.gg) the factor N has to be 
calculated for the finite lattices in use by taking into ac¬ 
count all details like choice of discretization, boundary 
conditions, the ratio c in ( |68] l etc. Note that in the SF 
there also is the indicated xq dependence and the obvi¬ 
ous standard choice is to take xq - T12 here. Depend¬ 
ing on the action in use N is calculated in closed form 
or numerically, involving a finite momentum sum, and 
the required values may be tabulated once and for all. 
The effort in any case is negligible. 

In figure we have a first look at the cutoff depen¬ 
dence of a step scaling function of a GF coupling a la 
El . Some details of this experiment are: Nf - Q,T - 
L,c - 0.3 with the actions S lw for the U sampling, S w 
in the flow equation ( [70| ) and a clover-discretized E. We 
see a nice extrapolation to the continuum limit but also 
a non-negligible slope and a clear break off from an ap¬ 
proximate asymptotic a^ behavior beyond ajL > 1/8 
(vertical dashed line). More detailed investigations of 
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Figure 6: Step scaling function of a gradient flow coupling for N{ = 0. 


the cutoff effects of GF couplings are underway and the 
improved flow equation ( [72| is presently simulated. 

6. Some details on simulation algorithms 

The whole physics programme reviewed in this ar¬ 
ticle depends on the availability of appreciable (paral¬ 
lel) computer resources and their optimal exploitation 
with efficient algorithms (see Il42ll for a review). While 
the resources mainly have to be provided by outside su¬ 
percomputer centers, algorithmic research and develop¬ 
ment is a major occupation of the physics collabora¬ 
tions. Usually useful ideas can only be advanced if the 
physics that one tries to extract is understood in detail. 
This section is therefore meant to reflect some of these 
efforts without being able to cover the whole held. 

6.1. Hybrid Monte Carlo 

Algorithms of the Hybrid Monte Carlo (HMC) class 
are at present the only known method that is able to cope 
with QCD at Nf > 0, although for Nf - 0 and for other 
bosonic lattice held theories much more efhcient algo¬ 
rithms exist. The problem is the nonlocal contribution 
of the fermion determinant in ( |20| l. 

The HMC Il43ll44l method is based on molecular dy¬ 
namics. There for each link U{x,ii) a conjugate mo¬ 
mentum n(x,ij) in the Lie algebra is introduced and a 
Hamiltonian is considered 

'H[n,U]^-Yj^(7THx,)i)) + S[U] (73) 

xp 

with some lattice action S[U] that plays the role of a 
potential here. We now consider an enlarged ensemble 
with a partition function given by the path integral 

Z= r (74) 


where the tt helds are integrated with the obvious mea¬ 
sure over their 8 dimensional real vector space for each 
link. It is clear that observables depending only on U 
assume the desired expectation values in this ensemble 
as the Gaussian n integrations simply factorize out. 

The HMC update procedure for the enlarged ensem¬ 
ble is now given by the following sequence of steps: 

• a complete field of independently Gaussian dis¬ 
tributed n(x,fj.) is drawn at random, 

• with the given configuration (tt, U) as initial values 
at f = 0 Hamilton’s equations are solved in molec¬ 
ular dynamics tim^f 


7T(x,fJ.) = -dx,pS[U] 

(75) 

U(x,fJ.) - 7T(x,fi)U(x,)l) 

(76) 


up to some trajectory length t - t. In practice 
this is done with some inexact discretized integra¬ 
tor with a hnite step size that has to be exactly time 
reversible, 

• only due to step size errors of the integrator the 
Hamiltonian between f = 0 and t - t will change 
HI ^ HI + AHl. The end configuration is taken as 
a successor of the initial one with the Metropolis 
acceptance probability min(l,exp(-A'7f)). In the 
case of rejection the old configuration remains un¬ 
changed. 

An important point here is that no zero step size limit is 
needed to prove detailed balance for the HMC. In prac¬ 
tice however the step size r/Astep has to be small enough 
to make the rejection rate small, around 10% for exam¬ 
ple in practice. 


This t has nothing to do with the flow time before. 









Rainer Sommer and UUi Wolff/Nuclear Physics B Proceedings Supplement 00 (2015) 7-|jO| 


14 


With two dynamical quark species for S an effective 
action like 

=5w[t/]-lndet(M2) (77) 


• for a pair of degenerate improved SW quarks there 
is no sign problem although the Dirac operator 
Dsw has complex eigenvalues. This is so since the 
determinant is real because of the relation 


should be taken that includes the fermion determinant. 
We have introduced here the usual lattice parameteriza¬ 
tion for the Dirac matrix (for one flavor) 

M = 2Ka(D -I- mo) (78) 

with the hopping parameter 


M* = ysMy^. 


(82) 


With the pseudofermion representation of the deter¬ 
minant the HMC Hamiltonian becomes 

'K[7r, = ^(n,n) + Sw[U]+Spf[U,(l>] (83) 


A. - , \ ^ J 

8 + lamo 

where D will be Dsw in most cases here. Its non¬ 
negative squared determinant corresponds to two de¬ 
generate quark species as in the Nf - 2 theory. Other 
cases complicate simulations further and will only be 
mentioned later. The advantage of HMC in this con¬ 
text is that in only the response of 5eff to in¬ 

finitesimal rather than finite moves in U, which enters 
into matrix elements of M, is required. The derivative 
of the Indet contribution in 5eff leads to the necessity 
to control arbitrary matrix elements of the non-sparse 
matrix M '. In four dimensional QCD on large lattices 
also this is impractical. Therefore the additional trick of 
introducing pseudofermions is needed that we discuss 
next. With them we shall only have to solve systems of 
linear equations with the Dirac matrix as coefficients for 
which highly efficient tools exist. 

6.2. Pseudofermions 

A pseudofermion is a complex field (j){x) which car¬ 
ries the same Dirac and color indices as one quark 
species. The squared determinant can now be repre¬ 
sented by a Gaussian path integral 

|det(M)|2 = jD[0]exp(-5pf[,^]) (80) 

with 

( 81 ) 

X 

A number of comments are in order; 

• D[0] means independent integration over the real 
and imaginary part of all components and some 
(later irrelevant) normalization, 

• due to the inversion of M the action of 4>{x) is 
highly nonlocal, 


with a scalar product notation for the n part and Wilson 
gluons as an example. Note that we have not introduced 
momenta conjugate to f. As they are Gaussian a cor¬ 
rectly distributed f configuration can be trivially drawn 
(for given U) by applying once the Dirac operator in 

(j) - Mr], (84) 

Here rj has the same indices (color, Dirac) as f and each 
component is an independent Gaussian random number 
with zero mean and unit variance. Now each HMC tra¬ 
jectory starts with a global choice of n and <p and then an 
evolution in molecular dynamics time t - 0 ,..., t and 
the accept/reject step. The force during the evolution 
stems from the gluon action and from 5pf. The deriva¬ 
tive of the gluon action for a given link is given 
by a sum of small Wilson loops, ‘staples’ for 5w, that 
require small computation time. The derivative of Spf 
derives from the variation of the discretized Dirac op¬ 
erator M ^ M + 6M as the U(x,fj.) entering its matrix 
elements change. It is given by 

X 

with 

Mx = (86) 

For each computation of the pseudofermion force we 
thus have to solve two linear systems involving the 
Dirac matrix. It is this step - solving the lattice Dirac 
equation with given right hand sides - which consumes 
by far the most time in QCD simulations with the HMC 
algorithm. 

Obviously it is worthwhile to optimize the linear 
equation solvers used with HMC as far as possible and 
a corresponding effort has been devoted to this issue 
by the lattice community. As the Dirac matrix in any 
of the discretizations discussed here is sparse, iterative 
Krylov space solvers are the method of choice. At first 
the simple and robust conjugate gradient method ll45l 
has been widely used. A large number of improvements 
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have been developed over the years with incremental 
small speed-up factors. A kind of breakthrough with 
regard to the slowing down for small quark masses has 
been achieved more recently with low mode deflation 
and with multigrid methods IH^ . 

Also the integrators used for Hamilton’s equations in¬ 
side of HMC allow for improvement. One example is 
the multiple time step size technique ll49l . It is based 
on the observation that the pseudofermion force, that 
is expensive to compute, is usually much smaller than 
the gluonic force. The former would therefore admit 
larger time steps leading to fewer evaluations. In ll49ll a 
modified leapfrog integrator is proposed that allows for 
different step sizes for the two components while still 
maintaining the required reversibility. Other possibili¬ 
ties are to minimize the higher order step size errors and 
a popular integrator in this respect has been proposed in 

Hi. 


• and IT, 1, 

• R can be represented by one or several pseud¬ 
ofermions. 


A well known solution is the rational approximation of 
degree m Ea 


R = 


m-1 




m\Ms -t- o)^ 

M'Im, + vl ■ 


(90) 


If one decides on a spectral interval of m]m, 

on which R approximates the inverse square root then 
there is the well defined Zolotarev algorithm to con¬ 
struct those C,{cok],{vk] that optimize the approxima¬ 
tion quality in a certain norm. In addition the shifts can 
be taken real, positive and ordered 


0 < Vo < Wo < Vl < ■ ■ ■ < w,„_i. 


(91) 


6.3. Hasenbusch preconditioning 

Hasenbusch lISTIl has proposed a simple factorization 
of the two flavor determinant 

I det(M)|2 = I det(M)|2 x | det(MM^')|2 (87) 

which has been further investigated in ||52| . Above, M 
contains a larger mass (k < k). One pseudofermion is 
introduced now for each factor with an action 

5pf = (M^Vi>^^Vi) + (^Af^V2,MM^V2)-(88) 

With some tuning of k the condition numbers of both 
factors can be lower than the one of the original M. As 
a consequence the forces are of smaller magnitude with 
smaller fluctuations and the step size can be increased 
by a factor two or so. It turned out that it was important 
to tune K properly 1531 and also more than two factors 
can be introduced. 

6.4. Simulation of nondegenerate quarks, Nf > 2 

All simulation techniques and improvements dis¬ 
cussed so far have exploited that the a fermion deter¬ 
minant was squared for two degenerate fermions. Its 
weight could not go negative and the factorization was 
useful for introducing a pseudofermion. To model QCD 
more precisely the inclusion of the strange (and ulti¬ 
mately also charmed) quarks becomes necessary. 

The strange quark implies an additional weight 
det(M,) to be included in the Boltzmann factor where 
the Dirac matrix M, includes the strange mass. Follow¬ 
ing H and H we start from the trivial factorization 

det(M,) = IT,det(/?^‘), W, = det(M,/?). (89) 

The goal is to construct R such that 


It is now not difficult to see that, using the rational fac¬ 
torization, det(/? ') can be represented by one or sev¬ 
eral pseudofermions. By attaching the subsets of the 
factors with similar shifts to several pseudofermions, 
even something similar to mass preconditioning can be 
achieved. 

In spite of a good rational approximation (at moder¬ 
ate degree m) we still have to take care of the correction 
factor Ws. Due to the lack of chiral symmetry the Wil¬ 
son type lattice Dirac operator has no rigorous spectral 
cutoff and may in principle even develop zero eigenval¬ 
ues for some gauge fields occurring in the path integral 
(Monte Carlo). In such cases IT, may not be close to one 
and therefore important. Such fluctuations are known to 
become extremely rare once the mass is not too small 
and additional light quarks are present. They are further 
suppressed for large volume. For the strange quark IT, 
should be monitored to have only moderate fluctuations 
around unity. In precisely this case it can be estimated 
stochastically by 

IT, = ^ 92 ) 

with Gaussian random fields 77. Its estimator, possibly 
averaging over several independent q fields, may then 
be included in observables as usual for reweighting. 

Clearly, stochastic strange quark reweighting can not 
account for fluctuations to negative det(M,). As indi¬ 
cated this is however not expected to happen for large 
enough mass and volume. In massless Nf - 3 simu¬ 
lations for renormalization purposes the finite volume 
SF boundary conditions have to be monitored via W, to 
sufficiently stabilize the spectrum. In the applications 
reviewed here this is indeed the case. 
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6.5. Topological freezing 

The efficiency of any Monte Carlo algorithm depends 
on the degree in how far it is able to produce at reason¬ 
able cost statistically independent lattice field configu¬ 
rations distributed according to the action of the field 
theory under study. The statistical error of any observ¬ 
able is given by 

a-HO)^v(Oy^. (93) 

In this relation viO) is the variance of O which is another 
expectation value. The number of generated estimates 
N is divided by the integrated autocorrelation time Ti^ip 
which summarizes the decorrelation power of the algo¬ 
rithms for a given observable. There also is a maximal 
autocorrelation time Texp that is related to the spectral 
gap in the transition probability (Markov) operator of 
an algorithm. More details on autocorrelations and er¬ 
ror analysis can be found in ll56l for example. 

Each smooth configuration in a continuum Yang 
Mills theory on a four dimensional torus carries an inte¬ 
ger valued topological quantum number Q. In the path 
integral all configurations are summed over and all Q 
contribute even if large \Q\ may have a small relative 
weight in some situations depending on the physical 
size of the torus and the observable. The famous self¬ 
dual instanton configurations, for example, are classi¬ 
cal fields with Q - +\. Quantization and in particu¬ 
lar the lattice regularization eliminates in the first place 
all naive arguments based on continuity and topology. 
It has been argued however Ell that topological wind¬ 
ing numbers do effectively imply a decomposition of the 
space of U(x,p) configurations into sectors that are sep¬ 
arated by barriers with large Euclidean action. The nu¬ 
merical problem now is that known Monte Carlo algo¬ 
rithms for QCD, including HMC, change configurations 
in small steps and the probability to change sectors by 
sequences of such moves gets very small when the con¬ 
tinuum limit is approached. As one lowers a this is first 
seen in steeply growing autocorrelation times. Once 
these are very large, they are hard to even discover in 
the statistical error analysis and there is a risk of un¬ 
noticed systematic errors. Clearly observables formed 
in terms of (a lattice approximation) of Q are optimal 
diagnostic tools to control this problem of topological 
freezing. 

The seriousness of the freezing problem was noticed 
during the physics programme reviewed here and led 
to quite some delay. A precise analysis of the problem 
was given in ll58l . Here a special strategy to sharpen the 
error analysis was proposed that allows for a reliable er¬ 
ror estimation in the range of lattice spacings down to 


a =» 0.05 fm where with our standard algorithms and 
actions freezing starts to become a problem. The idea 
of the method is to extract the longest autocorrelation 
time from measuring and E(t) at finite flow time and 
feed this information into the error analysis of other ob¬ 
servables that couple to topology only weekly. Slightly 
different strategies can be used for the running GE cou¬ 
pling with SE boundary conditions at intermediate vol¬ 
ume ll59l which however also suffers from freezing. 

Einally in l60l a fundamental solution to the problem 
was given. Open boundary conditions in Euclidean time 
- while keeping periodicity in the three spatial direc¬ 
tions - abolishes the quantization of topological charge 
and thus the barriers in field space. Charges can so 
to speak flow in and out of the boundary hyperplanes. 
In view of the seriousness of the freezing problem, the 
price of a somewhat more complicated data analysis in 
the large volume simulations due to the lost time trans¬ 
lation invariance seems reasonable to accept. In addi¬ 
tion Symanzik improvement has to be modified due to 
the presence of boundaries in a way that is somewhat 
familiar from the SE. 

Eor a recent review covering lattice QCD algorithms 
in more detail than possible here we refer to ISD. 


7. Running of coupling and masses 

In this section we review numerical results that have 
been achieved for the running of the non-perturbatively 
defined coupling and mass with the finite system size as 
renormalized scale in the Schrodinger functional. 


7.1. Step scaling function of the SF coupling at N{ - 2. 

The SE coupling has been defined in ( [47| and its evo¬ 
lution is analyzed by computing the lattice step scaling 
function (SSE) ( fTT) and by extrapolating it to cr. In 
the numerical investigations it has turned out that it is 
highly profitable to accelerate the continuum limit of E. 
We consider the perturbative expansion of the deviation 

Y,(u,alL) — cr(u) ^ 

-—- = 6 \{alL)u + 62 {alL)u^ + _(94) 

cr(u) 

Eor the simulations on which we report here the SE 
was implemented with the Wilson plaquette gluon ac¬ 
tion 5wsf with SW (‘clover’) improved quarks 5swsf- 
Eor this total action the two loop perturbation theory 
has been studied on the lattice llJTl and hence 61,62 are 
known for the relevant values of Lja. We use them to 
define a perturbatively two loop improved SSE 


lF\u,alL) 


'L{u, ajL) 

1 + 6 \{alL)u + 62 {alL)u^ 


(95) 
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Figure 7: Continuum extrapolations of the Schrodinger functional 
coupling with Nf = 2 dynamical quark flavors. Figure from (62). 

The very flat and well controlled continuum extrapola¬ 
tion of this quantity is shown in figure]^ Several other 
technical issues had to be mastered to produce these 
data. The two bare parameters go and the bare quark 
mass K for each Lja have to be tuned to the u val¬ 
ues of the series shown together with a vanishing quark 
mass. This can be achieved only to some limited preci¬ 
sion and small corrections have to be applied based on 
perturbative as well as numerical information. For the 
quark mass a particular definition m based on a PCAC 
relation ( [62| i is adopted and it is estimated that a 
tuning up to \mL\ < 0.05 sulfices for the attempted pre¬ 
cision. 

In figure we find an additional demonstration that 
the values of Lja - 6 ... 12 have very small discretiza¬ 
tion errors, at least after our 2-loop improvement of the 
observable and with the Wilson plaquette gauge action. 
One can therefore carry out a precise continuum limit 
with these rather small lattices. On the other hand the 
figure shows that one cannot take this for granted for 
any action. Care to take the continuum limit is the most 
important requirement for a trustable determination of 
the step scaling functions and ultimately also the A pa- 



Figure 8: A test of the continuum extrapolations with different actions 
for Nf = 0. The data from top (triangles) to bottom (open circles) ai‘e 
for the Iwasaki, the tree level Liischer Weisz and the Wilson gauge 
action. Both the boundary improvement of the action and the im¬ 
provement of the observables have been included. At present this is 
possible at the 2-loop level for the Wilson gauge action only, and at 
the 1-loop level in the two other cases. Figure from (SI based on data 
from I63II64I . 



Figure 9: The Af-dependence of the step scaling function of the 
Schrodinger functional coupling (m. Non-perturbative results are 
shown together with the two-loop curves. 


rameter. 

The continuum extrapolated step scaling function can 
finally be iterated to construct the non-perturbative run¬ 
ning coupling for a number of scale arguments in fig¬ 
ure We will see that changing the number of quark 
flavors from Af = 0 to Af = 2 does not induce any 
qualitative changes. The connection from low to high 
energies is rather smooth and perturbation theory can 
be trusted in the Schrodinger functional schemes rather 
precisely at energies >50 GeV or larger. For the pur¬ 
pose of the comparison we show side by side Af = 0 re¬ 
sults from where the methods were developed and 
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Figure 10: Running coupling for A^f = 2 compared to A^f = 0. Figure from (29) based on results of I66II67|[^I64II62I . 


Nf - 2 results. 

With the continuum step scaling function cr{u) under 
our control for a range u e (0, Umax) we can connect the 
associated scales Asp and L^ax implicitly defined by the 
condition 

g^iLmax) = Mmax = 4.484. (96) 

The result, including a complete error analysis, is ll69]l 

TmaxA^^i^ = 0.264(15). (97) 

In terms of the semi-phenomenological scale ro ~ 
0.49fm the estimate L„ax ~ 0.39fm can be given. This 
scale is defined llTOll in terms of the force F between 
static quarks by solving 

rlF(ro) = 1.65. (98) 

It is loosely connected to quarkonia models and has in 
addition been related to other observables in previous 
lattice simulations. In the next section more direct con¬ 
nections of Lmax to phenomenology will be cited. 


7.2. Verification of asymptotic freedom 

Asymptotic freedom is normally taken for granted 
for QCD with any 0 < Af < 6 and it is certainly a 
self-consistent property in perturbation theory. Once 
we have some control beyond perturbation theory we 
should however remember this situation and analyze our 
data under this aspect. 

From the continuum limit SSF we may form the com¬ 
bination 


bf{u) = 


cr{u) - u 
2 log 2 


(99) 


It is expected to extrapolate to the perturbative coeffi¬ 
cient bo at small u, with an asymptotic approach that 
is linear in u. In figure where we include additional 
data for Af = 0,2,3,4, we do see the expected behav¬ 
ior. This provides a clear non-perturbative confirmation 
of asymptotic freedom including the universal perturba¬ 
tive Af dependence at small coupling. 

In perturbation theory asymptotic freedom is lost be¬ 
yond Af = 16 because bo changes its sign. We just re¬ 
mark in passing that we see a considerable effort to de¬ 
termine SSFs for Af = 8 ... 12 where one expects an al¬ 
most vanishing jS-function and approximate conformal 
invariance that may be of phenomenological interest for 
models of the technicolor variety. Many of these inves¬ 
tigations use the methods described in this article, but in 
a much more difficult situation. 


7.5. Running quark mass 


We now mention also some result for the running of 
the quark mass which is actually intertwined with the 
evolution of the coupling constant. The scale depen¬ 
dence of the running mass m(L) in (|6^ derives form the 


L dependence of Zp in (60 1 . The renormalization group 
equation for the mass reads 


p— = T(g)m, L ^ (100) 

d/j. 

with T{g) - -%g^ HAnf + ©(g"*) in perturbation theory. 
Note that this scale evolution is coupled with the one of 
giif in Q- To have non-perturbative control also here. 
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Figure 11: Continuum extrapolations of the step scaling function Ep in the quenched approximation (left) (13 and with Nf = 2 dynamical quark 
flavors ED. In the right graph the coupling u ranges from u = 0.979 to w = 3.33. 


we define another step scaling function 
Zp(go, 2L/a) 


Sp(M, ajL) - 


Zp{go,Lla) 


( 101 ) 


u=f{L) 


Also the renormalization group equation for the mass 
can be converted to an equivalent integral equation 


T(g) 


( 102 ) 


exp 


/ 


dg 


pig) 


do 

bog 


where M is the scale independent RGI mass, that is the 
analogue to the A parameter in the case of the coupling. 
The perturbative approximations of t and p may be used 
at very high energy (small L in the SF) to compute M, A 
from m{L) and giL). After extrapolating to the contin¬ 
uum limit 


crp(M) = lim l,p(u,alL), (103) 

a/L^O 

see figure El the continuum running mass was con¬ 
structed. It is shown in figure [T2] 


7.4. Optimized strategy 

In figure we have already shown some results for 
the running of the SF coupling with Nf - 3. The goal in 
the present stage of the project is to improve the preci¬ 
sion at the same time as having this more realistic num¬ 
ber of flavors. As already mentioned at the beginning 
of section]^ the combination of SF and GF coupling is 
promising for the precision issue. In figure [T^ we show 



Figure 13: Sketch (no real data) of the overall strategy being employed 
forWf = 3 (m. 


a sketch of the overall strategy to use their complemen¬ 
tarity and combine the two couplings. We note however, 
that here we do not yet see real data! 

The presently quite different size of the discretiza¬ 
tion errors in the two observables are seen by compar¬ 
ing figure and figure We hope that the improved 
flow equation ( |7^ and the other optimizations reviewed 
in that subsection will still considerably accelerate the 
continuum limit of the GF coupling. 

Increasing the precision also requires a more precise 
tuning to the massless theory than before. This is at the 
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m(M)/M 



Figure 12: Scale dependence of the quark mass m in the quenched approximation tleft) l67l and with A^f = 2 dynamical quark flavors ITll . For 
A^f = 0 the dotted, dashed and solid curves are obtained from eqs. and GD using the 1/2-, 2/2- and 2/3-loop expressions for the r- and 
yS-functions respectively. 
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Figure 14: Determination of the critical lines for various values of 
Lja. Note that a big part, the 2-loop expression, is subtracted from 
the data. Figure from EH. 


same time desirable for other projects such as HQET, 
where the precise knowledge of the relation between 
g and L is used to perform simulations in a volume of 
fixed physical size. For SW improved Wilson fermions 
that we use this means that the critical mass in lattice 
units amdgo) has to be determined to a many digit pre¬ 
cision in the relevant range of couplings. For the Wilson 
gauge action such a result is given in figure For fu¬ 
ture simulations these data are represented by a smooth 
interpolating fit formula shown as lines in figure [T4| 


8. Large volume simulations 

In this section we review the completed results for 
Nf - 2 and some data of the ongoing Nf - simula¬ 
tions. Quenched computations from the nineties were 
an extremely useful preparation but will be skipped 
here. 


8.1. Algorithmic issues 

Here simulations are described that have been run for 
Nf -2 improved Wilson quarks within the Coordinated 
Lattice Simulation (CLS) consortium. 

The algorithmic framework is mostly as discussed in 
section where it has become clear that the frequent 
inversion of the lattice Dirac matrix is the dominant nu¬ 
merical task. It is in particular required to evaluate the 
fermionic driving force for the molecular dynamics evo¬ 
lution in HMC. 

Before we have switched to the more efficient Hasen- 
busch preconditioning (see subsection |6.3| l, the domain 
decomposed DD-HMC algorithm ll72l had been used 
initially. This at first very promising variant was ul¬ 
timately dismissed because of inferior autocorrelations 
and critical slowing down, but we here still look at some 
data produced with DD-HMC for comparison and for 
historical reasons. In figure 15 we see several histo¬ 
ries in molecular dynamics time, where the left plots 
refer to DD-HMC simulations and the right ones (MP- 
HMC) to Hasenbusch mass preconditioning. Compo¬ 
nents of the fermion forces from two separate pseud¬ 
ofermions are shown and the Hamiltonian violation A'TT 
that enters into the accept step. The spike-like fluctua¬ 
tions that one sees point to algorithmic problems due to 
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Figure 15: Histories of the energy violation AH, as well as maximum and average forces F 2 and Fi, for each force update, plotted as a function 
of the trajectory number. Values con'esponding to the DD-HMC algorithm are shown to the left and the integration step-sizes for the two forces 
relate as At 2 : At\ = 1:6. The values for MP-HMC are shown in the two right panels and the corresponding ratio of the integration steps is 
Afi : At 2 = 1:10. The lattice size is 48 X 24^ and /tsea = 0.13625. 

nearly singular Dirac matrices. With them present, one 
is restricted to small step sizes and finds large iteration 
numbers for the Dirac solvers which both damage the 
efficiency. The plots demonstrate the smoother running 
of the MP-HMC after tuning its parameters in a reason¬ 
ably close to optimal way GS. After detailed studies 
MP-HMC was adopted as the method of choice. 

Topological freezing (see subsection |6.5| l was the 
other problem that significantly held back the comple¬ 
tion of large volume Nf - 2 simulations. In the history 
of the topological charge in figure one sees that the 
simulations are slowing down dramatically as the lat¬ 
tice spacing is reduced. The observed structures extend 
over a non-negligible fraction of the length of a typi¬ 
cal run. On the other hand, studying the effect of the 
slow modes of the Monte Carlo algorithm on typical 
hadronic correlation functions, it was found that these 
receive only suppressed contributions to their autocor¬ 
relation functions. Therefore, for a ^ 0.05 fm the stan¬ 
dard etTor computation could still be adapted to the sit¬ 
uation 1581 ED. Conservative error estimates are still 


possible with MC histories which are of the order of 
20-100 times the slowest relaxation time of the system, 

Texp- 

8.2. General considerations for scale setting 

In the previous section the energy scales from the per¬ 
turbative end down to the implicitly defined hadronic 
scale Lmax has been covered for coupling and quark 
mass running. It remains to connect Lmax to a mea¬ 
surable quantity to bring in MeV units which we call 
‘scale setting’. After the general discussion in section]^ 
in principle, if lattice QCD had all dimensionless scale 
ratios right, any observable would be suitable and the 
mass of the stable proton would appear particularly nat¬ 
ural. In reality, however, many more mundane practical 
considerations are essential GD to not give away too 
much precision in this step. 

An obvious demand is that the quantity used for scale 
setting should be computable on the lattice with good 
statistical precision (for a given computational effort) 
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a=0.08fm a=0.06fm a=0.04fm 





Figure 16: Histories of an estimate of the topological charge on runs with A^f = 2(74). 



Figure 17: Elfective masses for nip Gl], mn ES), V{== ro), y(« n) (69) and (77) on CLS ensemble N6 (see (77)). All elfective “masses” have 
been scaled such that the errors in the graph reflect directly the errors of the determined scales. They have been shifted vertically. Figure from (H. 


and that it should have small systematic errors like cut¬ 
off and hnite volume effects. The latter would otherwise 
contaminate all quantities cited in MeV even if in lattice 
units they themselves are well controlled. Another point 
is that it is desirable to use a scale that is not very sensi¬ 
tive to the precise quark content and to the tuning of the 
masses of included sea quarks to their physical values. 
This is obviously only possible to a limited precision 
and in addition, for algorithmic reasons, the up/down 
quarks in simulations are most of the time heavier than 
in Nature and a chiral extrapolation is invoked. This 
step is theory guided by chiral perturbation theory, but 
it is clearly good to keep the general scale as indepen¬ 
dent of this as possible. In addition the effective theory 
with Af < 6 is an approximation only and scale setting 
should be carried out in the sector of the theory that is 
robust against these presently still unavoidable small er¬ 


rors. 

With this said, estimates for a number of quantities 
that are considered by the community are displayed in 
hgurefT?! The plot includes the effective masses of the 
proton and the Q particle made of three strange quarks 
(still quenched for Nf = 2). We see the relatively large 
errors combined with short plateaus, where the situa¬ 
tion for the Q. is somewhat superior. The static quark 
potential V{rQ) is closely related to ro in ( [98| ) and y(ri) 
is a similar quantity where 1.65 is replaced by 1.00. The 
lowest line with a very convincing plateau with small er¬ 
rors represents an effective matrix element correspond¬ 
ing to the decay constant f„. A further bonus of such 
quantities is that they have been confirmed to be only 
weakly coupled to the slow modes responsible for topo¬ 
logical freezing. Therefore a reliable error estimation 
is still possible down to the smallest lattice spacings 































Rainer Sommer and UUi Wolff/Nuclear Physics B Proceedings Supplement 00 (2015) l-\30\ 


23 


a ss 0.05 fm entering here. A small draw back asso¬ 
ciated with decay constants has to be mentioned too: to 
cite a physical value for them based on physical decays, 
separate experimental input for the relevant CKM ma¬ 
trix elements has to be used. 

What has been said about f„ is also true for the Kaon 
decay constant /k. It is our preferred scale setting quan¬ 
tity in this study, because of two more bonuses. The chi¬ 
ral extrapolation from larger pion masses in the range of 
around 500 MeV to 270 MeV down to the physical value 
is simpler for /k. In this case the partially quenched 
variant of chiral perturbation theory (pqChPT) is used. 
A purely technical point to prefer the K sector in our 
Nf - 2 simulations is that the strange quark mass can 
here be varied after the tTin without having to gener¬ 
ate new configurations. Only a few extra inversions for 
the strange quark propagator are needed which is a rel¬ 
atively small effort. 


8.3. Scale setting for Nf — 2 

We start with a compilation of simulation data in ta¬ 
bles]^ an d[T] The non-integer values for Lmax/a derive 
from interpolations, see 1 ^ for more details. 

The main difficulty and source of a systematic error is 
the extrapolation to the proper quark masses, the “phys¬ 
ical point”. Once we decide to set the scale through /k, 
this point is naturally defined by 


Rk^R' 


,phys 


R.^R' 


,phys 


where 


2 2 
- TP > - 72 ’ 

-'K -'K 


(104) 


(105) 


and are the values of these ratios in Nature. 

In an attempt to minimize uncertainties, we take the 
physical masses and decay constants to be the ones in 
the isospin symmetric limit with QED effects removed 
as discussed in II 8 OI . We use 


'«;r,phys = 134.8MeV,mx.phys =494.2MeV. (106) 

One can then, for each lattice spacing, carry out 
two different strategies for extrapolating to the physical 
point. 

In strategy 1 one keeps 

R^="^^R^^y\ (107) 

7k 

fixed, as one varies the light (dynamical) quark mass. 
The condition defines a curve in the quark mass plane 
spanned by (m^ -)nia and m^. It is a very interesting one 


because along that curve tok is constant up to terms of 
order mj. In the ChPT expansion the mass-dependence 
is then small and in particular the coefficient of the chi¬ 
ral log, TO^logmJ, is small. One expects that /k can 
be extrapolated rather easily to the physical point along 
this curve. 

The order corrections are known in terms of one 
low energy constant, ua,. One just has to implement our 
condition eq. (107 1 , which expresses otk in terms of m„, 
in the formulae of EU- 

The predicted form is 


/k 

- /K.phys [1 + TK(yi,yK) 

(108) 


+{aA - 7 ) (fi - yif + 0 (y^), 

i'K(fi.yK) 

- i'K(yi>fK) - i'K(y;r>yK) . 

(109) 

i'KCVl.fK) 



log(2yK/yi - !)■ 

( 110 ) 

The variables 






yn 


JK 


’ 


u 

;r,phys 




K,phys 


K,phys 




K,phys 


0.00958, 

0.12875, 


are proportional to (averages of) quark masses up to 
quadratic terms. Because of eq. ( 107| i, we have ^3 s 
m^/[ 87 r^/^] = 2 yK - yi + O(y^) and ^3 does not appear 
in eq. (108 1 . 

Another option is strategy 2, where one keeps the 
(PCAC) strange quark mass fixed and uses the expan¬ 
sion in just the up quark mass, i. e. SU(2) chiral pertur¬ 
bation theory adopted to this situation. 

Figure [TS] shows extrapolations with both strategies, 
which converge well at the physical point. 

It remains to perform a continuum extrapolation of 
the dimensionless combination /Kfmax with the help of 
interpolations of the integer Lmax/a as a function of the 
bare coupling g^. Little discretisation errors are seen in 
figure [T^which leads to the continuum limit 


/KL,nax=0.315(8)(2). (Ill) 


As in the previous section, the final results come from 
strategy 1 for the chiral extrapolation of Fk. Strategy 2 
is used to estimate the systematic uncertainty in the sec¬ 
ond parenthesis; it is small compared to the statistical 
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L/fl 

P 

^sea 

am 

fiD 


^max/^ 

4 

5.2000 

0.134700 

-0.03745(41) 

3.730(11) 



4 

5.2000 

0.133780 

-0.00086(35) 

3.797(11) 



4 

5.2000 

- 

^ 0 

3.798(11) 

-0.686 

5.11(3)(13) 

6 

5.2000 

0.135600 

-0.01322(26) 

4.810(32) 



6 

5.2000 

0.135200 

+0.00289(24) 

4.984(33) 



6 

5.2000 

- 

^ 0 

4.954(33) 

+0.470 

5.33(4)(11) 

6 

5.2638 

0.135673 

+0.00012(19) 

4.550(25) 

+0.066 

5.89(4)(2) 

8 

5.4689 

0.136575 

+0.00046(11) 

4.526(32) 

+0.042 

7.91(7)(1) 

10 

5.6190 

0.136700 

+0.00038(8) 

4.531(51) 

+0.037 

9.87(14)(2) 

12 

5.7580 

0.136623 

+0.00067(7) 

4.501(91) 

+0.017 

11.94(31X1) 

16 

5.9631 

0.136422 

-0.00096(4) 

4.40(10) 

+0.084 

16.40(50)(6) 


Tsblc 1. Values of flftcr coiT6ctiri§ the simul^tcci vslucs Bja to ni3.tch the g from (69) The data at the two lai'gest yS-values are from 

ED. The second error on the final result is the systematic one. 


P 

^max/^ 

imax/K 

rolL^ax 

^Strange / /k 

5.2 

5.367(82) 

0.318(6)(3) 

1.155(22) 

0.530(12)(6) 

5.3 

6.195(51) 

0.320(5)(4) 

1.169(15) 

0.577(11)(7) 

5.5 

8.280(80) 

0.316(4)(2) 

1.213(17) 

0.617(11)(5) 

cont. 


0.315(8)(2) 

1.252(33) 

0.678(12)(5) 


Table 2: Values of Lm^xlo, Lmax/K. ^o/f-max and mstiange//K together with the values extrapolated to the continuum limit. The running mass in the 
Schrodinger Functional scheme mstrange is given at the renormalization scale Z^^ax- 



Figure 18: Physical point extrapolation of the kaon decay constant in 
lattice units. Open symbols and dashed lines correspond to strategy 
1, whereas filled symbols and dash-dotted lines represent strategy 2. 
Only data below yi =0.1 enter the extrapolation. Figure from 


errors. We then quote 

A^i^//K = 0.84(6). (112) 


Now, as a result of our analysis, the error is dominated 
by the error on AL^ax- We translate to the MS scheme 
using = 2.382035(3)A®|I82] |83l as well as to 
physical units 


A:^ = 310(20)MeV, 


where 


/K,phys - 155 MeV 


(114) 


enters. 

As discussed previously with not all flavors of QCD 
treated dynamically there is a small ambiguity in the 
translation to MeV. We therefore also give the result 

ro A® = 0.789(52), (115) 

based on 


ro/imax = 1.252(33) 


(116) 


in complete analogy to eq. ([TTT]i. Our result 
an unambiguous non-perturbative property 
flavor theory, sometimes called QCD-lite. 


eq. ( |115| l is 
of the two- 


8.3.1. Strange quark mass 

The RGI mass Ms is given in terms of the bare PCAC 

mass ^strange 6y 

M 

lUg — ^strange (T) (117) 

M Za 

“ m(L)Zp(L)'”“’ 


(113) 


(118) 
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Figure 19: Continuum extrapolation of Lj = L^ax in units of /k- 
Even though the data shows no cut-off effects, a linear extrapolation 
is used to account for uncertainties from O(fl^) effects hidden by the 
eiTors. The two strategies for the chiral extrapolation of Fk = a/k 
agree well within statistics. 


where mstrange is the PCAC strange quark mass of strat¬ 
egy 2. The continuum value of the universal first factor 
M/m has been computed in Ref. ED and discussed in 
section |73] 

Expressing Mg in units of /k we eliminate Za and get 


strange 

/x.phys 


M 

m{L) 


^^strange 

X —:-X 

rbare 

^ K,phys 


(119) 


Zp(L) 


[1 + (^A ^p)^^strange] i 


with = E'K,phys/ZA- The second factor is 0(a) 

improved, if we neglect a tiny correction proportional to 
the sea quark mass, (Ija - bp)amsea- Note that Da, Dp = 
O(gQ) are loop-suppressed and is very small. 

Our final result is 


mstrange//K = 0.678(12)(5) , (120) 

Mg//K = 0.887(19)(7), (121) 

Mg = 138(3)(l)MeV, (122) 


where we use M/m = 1.308(16) at the scale L^ax- For 
reference, we also give the numbers in the MS scheme. 
This conversion is the only part of the computation in 
which we need to take recourse to perturbation theory, 
known in this case to four loops Il84lf85ll8^[87l . which 
differs from the two- and three-loop result by only a 
small amount. We use the same method as described 
in ED, but with the new value of Aj^ which leads us 
to m^(2 GeV)/M = 0.740(12) and 


mf®(2GeV) 


Mg m“^(2GeV) 

7k M 
102(3)(l)MeV . 


/x.phys 


(123) 

(124) 


8.4. Scale setting for Nf - 2 + \ 

After the progress in simulation algorithms and the 
understanding of how to get around the topological 
freezing, CLS has started large-scale QCD simulations 
with a strange quark in addition to degenerate up and 
down quarks. The action 5LW;SWsf is used with non- 
perturbative csw ED- The simulations started just about 
15 years ago, but have already reached a similar cov¬ 
erage of lattice spacings and pion masses as the Af = 
2 simulations carried out before. A summary of the 
presently available ensembles is found in table from 
|[8^ where the details of the simulations are described. 

With a dynamical strange quark, the choice of the 
curve in the quark mass plane that one chooses to ap¬ 
proach the physical point (cf. Sect. |8.3| l is much more 
important since each choice of the strange quark mass 
means a new simulation. Various considerations enter 
into a choice of a trajectory. Strategy 2 of Sect. |8.3| 
seems a natural one, but one has to keep the renormal¬ 
ized strange mass fixed and in the non-degenerate case, 
the renormalization of the quark mass contains a mix¬ 
ture of the flavor singlet and the flavor non-singlet mass 
terms, which renormalize differently. Similarly the 0{a) 
terms become more complicated ll89ll . As a result of this 
it is technically simpler to keep the trace of the quark 
mass matrix constant as one changes the light quark 
mass EOl. This condition is similar but not the same 
as strategy 1. The most important simplification is that 
up to a (supposedly small) 0(a) term, the condition of 
a fixed trace of the renormalized mass matrix is equiva¬ 
lent to 

■sp 1 

> —= const. (125) 


It can thus be followed without any tuning errors. A 
non-trivial point is of course to choose the right value of 
the trace, the one which leads to a trajectory through the 
physical point. Slightly wrong choices and subsequent 
corrections are unavoidable. We skip this issue here. 


The scale setting will proceed in analogy to Sect. 8.3 


through the decay constants With open bound¬ 

ary conditions translation invariance in time is lost. 
Boundary-effects of correlation functions exist and have 
to be taken into account. The theoretical analysis of 
these effects is the same as with SF boundary conditions 
in a large volume ll9T1l . Numerical aspects are presently 
being studied in detail 1921IMI . A short summary is that 
the boundary conditions do not present an obstacle for 
the extraction of the hadronic matrix elements such as 
decay constants. In fact for some cases they may be ad¬ 
vantageous compared to the conventional torus, where 
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id 


iVs 


Ku 

Ks 

m„[MeV] 

raxiMeV] 

m„L 

BIOS 

3.40 

32 

64 

0.136970 

0.13634079 

280 

460 

3.9 

HlOl 

3.40 

32 

96 

0.13675962 

0.13675962 

420 

420 

5.8 

H102 

3.40 

32 

96 

0.136865 

0.136549339 

350 

440 

4.9 

H105 

3.40 

32 

96 

0.136970 

0.13634079 

280 

460 

3.9 

ClOl 

3.40 

48 

96 

0.137030 

0.136222041 

220 

470 

4.7 

DlOO 

3.40 

64 

128 

0.137090 

0.136103607 

130 

480 

3.7 

H200 

3.55 

32 

96 

0.137000 

0.137000 

420 

420 

4.4 

N200 

3.55 

48 

128 

0.137140 

0.13672086 

280 

460 

4.4 

D200 

3.55 

64 

128 

0.137200 

0.136601748 

200 

480 

4.2 

N300 

3.70 

48 

128 

0.137000 

0.137000 

420 

420 

5.1 

N301 

3.70 

48 

128 

0.137005 

0.137005 

410 

410 

4.9 

J303 

3.70 

64 

192 

0.137123 

0.1367546608 

260 

470 

4.1 


Table 3: List of present CLS ensembles with up, down and strange sea quarks and the action of ED . The numbers for m^^ and mx are rounded 
and use = 0.4144 fm. The lattice spacings are roughly a = 0.086 fm, 0.064 fm and 0.05 fm foryS = 3.4, 3.55 and 3.7, respectively. Table from 

ED- 


0.11 
0.1 
0.09 

4 0.08 

^ 0.07 

0.06 
0.05 
0.04 

0 10 20 30 40 50 60 70 80 

xo/a 

Figure 20: Effective pion mass of the D200 lattice (see table[^. Graph 
from from Ell. 

particles may propagate around the periodic time. With 
open boundary conditions such effects are avoided. 

We show an example of an open boundary condi¬ 
tion correlation function in figure]^ from ll88l . Apart 
from the boundary effects the effective mass plot ex- 
hibts significant wiggles at large time, but these are of 
a purely statistical nature. Correlations between neigh¬ 
boring time-slices are very strong, but over larger dis¬ 
tances also anticorrelations are present and in the end 
the htted two-state curve is statistically compatible with 
the data points. 

For the extraction of pseudo-scalar decay constants, 
just as for their mass, one has various possibilities to 
combine correlation functions. There numerical study 
revealed that the details do not matter too much and 


given a certain number of decorrelated conhgurations 
the statistical precision of these quantities is comparable 
to the one with periodic boundary condition ensembles 

iaa. 

In summary, the 2 h- 1 simulations advance very fast 
and the scale setting is expected to be available rather 
soon. 

9. Outlook 

Let us first go back to hgure Here the present 
knowledge of the running of the SF coupling is summa¬ 
rized in an easily accessible way. For Af=3 the range 
of couplings is presently restricted. However, within 
this range the precision is already much better than for 
smaller flavor numbers. Still, the ALPHA collaboration 
is reducing the errors further. In parallel, work is now 
progressing on extending the range towards small ener¬ 
gies. As explained earlier, in the lower energy region it 
becomes advantageous to use the flow coupling. Now 
that the reduction of cutoff effects of flow observables 
by Symanzik improvement is understood BOll . we are 
ready to use it. The first preparation, the determination 
of the critical lines for different Lja in the necessary 
region of larger go as compared to figure is already 
far advanced. The step scaling functions of gcF will be 
evaluated soon. 

As we have reported in the previous section, also the 
necessary large volume simulations are far advanced. 
Hence, we foresee to soon present the three flavor A_ 
with a precision that is at least comparable to the present 
FLAG average 0. At this point, a perturbative rela- 
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tion of to can be used and one obtains an es- 

MS MS 

timate of a^^(Mz) with a precision of the non-lattice 
PDG average ||93]| or better. The additional uncertainty 
introduced by this step can roughly be divided into two 
pieces (for a precise description of the decoupling of 
heavy quarks on the non-perturbative level we refer to 
0). There are power-suppressed 0((Mchann/A|^) 
terms due to the neglect of higher dimensional opera¬ 
tors in the effective Lagrangian when we treat the low 
energy theory by just three flavor QCD. Such effects are 
entirely due to charm quark loops and are therefore sup¬ 
pressed by a factor of the number of colors in the large N 
expansion on top of a perturbative suppression, i. e. by 
a/N. In a (quite realistic, we would say) model it has re¬ 
cently been shown that indeed these power-suppressed 
terms are very small. We can safely neglect them within 
the envisaged accuracy. What remains is the matching 
of QCD with 3 flavors to QCD with 4 flavors. In terms 
of the A-parameters this is the relation 

A^ = P 3 4 (Mcharm/A—) A'^ , (126) 

which has a perturbative expansion in terms of the cou¬ 
pling at the scale jj. — Mcharm- In the MS scheme 
the relation is known to four loops and the resulting per¬ 
turbative uncertainty looks very small ll9^ . see 0 for 
the discussion of ^’(Mcharm/Aj;^). 

Nevertheless, perturbation theory at the scale 
ji - Mcharm is Worrying per se. Therefore the ALPHA 
collaboration also foresees a further step to carry out 
an adapted version of the full programme with four dy¬ 
namical quark flavors. A first step is to bring Symanzik 
0 (a) improvement under control with a heavy charm 
quark. We plan to carry out the steps at low and inter¬ 
mediate energy in a massive renormalization scheme 
mi. Concerning improvement, this scheme does not 
need to be defined exactly with the charm mass at its 
physical value, but it is sufficient to have it fixed and 
close to it such that in an expansion in amcharm - 
one can safely neglect higher order terms. This is not 
an easy undertaking, but first steps are promising ll95l . 
We can hence foresee to have in the near future a full 
four flavor non-perturbative determination of the A 
parameter. 
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